#define CATCH_CONFIG_MAIN
#include "catch.hpp"
#include "coherentElastic/coherentElastic_util/sigcoh.h"
#include "generalTools/testing.h"
#include <algorithm>


TEST_CASE( "full sigcoh" ){
  GIVEN( "Specify vec1 and vec2" ){
    double eFirstBragg = 4.555814e-4,
           scon = 321038.08426486229,
           temp = 296.0,
           recon = 5.1803120897E-20,
           emax = 0.01,
           eNext;

    std::vector<double> 
    vec1 { 8.794479E15, 8.794479E15, 3.517792E16, 3.517792E16, 7.915031E16, 
    7.915031E16, 8.717302E16, 8.717302E16, 8.717302E16, 9.596750E16, 9.596750E16, 
    9.596750E16, 1.223509E17, 1.223509E17, 1.223509E17, 1.407117E17, 1.407117E17, 
    1.663233E17, 1.663233E17, 1.663233E17, 1.930386E17 },
    vec2 { 0.000000E00, 0.000000E00, 2.0984634382E-8, 2.098463E-8, 0.000000E00, 
    0.000000E00, 1.626950E-9, 1.626950E-9, 1.626950E-9, 9.266134E-9, 9.266134E-9, 
    9.266134E-9, 2.702515E-9, 2.702515E-9, 2.702515E-9, 9.995415E-9, 9.995415E-9, 
    6.814547E-9, 6.814547E-9, 6.814547E-9, 6.814547E-9 };

    std::vector<double> epVec { 2*eFirstBragg, 5*eFirstBragg, 10*eFirstBragg};
    std::vector<std::vector<double>> sVariousEp{
      { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }, 
      { 5.914962,-3.548978, 0.236601, 2.129384,-2.413305, 0.9028636 },
      { 3.301423, 0.253589,-0.975299,-1.136607, 0.972065, 0.650655 } };
    std::vector<double> nextBraggEdge { 1.822325E-3, 4.100232E-3, 4.971415E-3 };
    
    WHEN( "Some number of legendre order are requested" ){
      for ( int nl = 1; nl <= 6; ++nl ){
        std::vector<double> s(nl,0.0), sThisEp(nl);
        THEN( "That number of s vector entries are returned, as well as the "
              "energy of the next bragg peak" ){
          for (size_t i = 0; i < epVec.size(); ++i){
            copy(sVariousEp[i].begin(), sVariousEp[i].begin()+nl,sThisEp.begin());
            eNext = computeCrossSections( epVec[i], vec1, vec2, emax, scon, recon,
                                      s );
            REQUIRE( eNext == Approx(nextBraggEdge[i]).epsilon(1e-6) );
            REQUIRE( ranges::equal( s, sThisEp, equal_1e5 ) );
          }
        } // THEN
      } 
    } // WHEN
  } // GIVEN
  
  GIVEN( "that vec1 and vec2 are generated using the bragg function" ){
    std::vector<double> vec1 (5000,0.0), vec2 (5000,0.0);
    int lat = 1, numAtoms = 1;
    double temp = 296.0, 
           emax = 5.0,
           recon = 5.1803120897E-20,
           eNext;
    double scon = prepareBraggEdges(lat, temp, emax, numAtoms, vec1, vec2);

    using std::copy;

    std::vector<double> epVec {5e-3,1e-2,5e-2,1e-1,5e-1,1,5};
    std::vector<double> nextBraggEdge { 0.0063381, 0.0113895, 0.0500739, 
      0.1014105, 0.5212158, 1.0063780, 5.0000010 };
    std::vector<std::vector<double>> sVariousEp {
      { 4.793003,-1.286693,  0.822574,-2.6573778,  1.8506794,-0.4372194 },
      { 3.954885, 0.0384121,-0.414184, 0.0202615, -0.5485786,-0.2275314 },
      { 4.271526, 0.147108, -0.104281,-0.177182,  -0.1513574, 0.2949411 },
      { 3.556954, 0.470383, -0.059366,-0.086328,  -0.0229501,-0.0795200 },
      { 1.209989, 0.689607,  0.263473, 0.062955,   0.0016168,-0.0175757 },
      { 0.610934, 0.473528,  0.289359, 0.141313,   0.0540541, 0.0132590 },
      { 0.122203, 0.116700,  0.106422, 0.092671,   0.0770404, 0.0611182 } };

    WHEN( "Some number of legendre order are requested" ){
      for ( int nl = 1; nl <= 6; ++nl ){
        THEN( "That number of s vector entries are returned, as well as the "
              "energy of the next bragg peak" ){
          std::vector<double> s(nl,0.0), sThisEp(nl);
          for (size_t i = 0; i < epVec.size(); ++i){
            copy(sVariousEp[i].begin(), sVariousEp[i].begin()+nl,sThisEp.begin());
            eNext = computeCrossSections( epVec[i], vec1, vec2, emax, scon, recon, 
                                      s );
            REQUIRE( eNext == Approx(nextBraggEdge[i]).epsilon(1e-6) );
            REQUIRE( ranges::equal( s, sThisEp, equal_1e5 ) );
          }
        } // THEN
      } 
    } // WHEN
  } // GIVEN


  GIVEN( "various materials are considered" ){
    std::vector<double> vec1 (5000,0.0), vec2 (5000,0.0);
    WHEN( "choose beryllium metal" ){
      int lat = 2, numAtoms = 2;
      double temp = 450.0, 
             emax = 6.0,
             recon = 5.1803120897E-20,
             eNext;
      double scon = prepareBraggEdges(lat, temp, emax, numAtoms, vec1, vec2);
  
      using std::copy;
  
      std::vector<double> epVec {5e-3,1e-2,5e-2,1e-1,5e-1,1,5};
      std::vector<double> nextBraggEdge { 0.00521980, 0.01159117, 0.05334958, 
        0.10194200, 0.50731040, 1.04618600, 5.05927600 };
      std::vector<std::vector<double>> sVariousEp {
        { 0.0000000, 0.0000000, 0.00000000, 0.0000000, 0.0000000, 0.0000000 },
        { 2.9931375,-0.9115081,-1.02714765, 1.0923217, 0.1193862,-0.8561407 },
        { 2.4814615, 0.0644388,-0.03095040,-0.2134906,-0.1475377,-0.0424251 },
        { 1.7162921, 0.4093202,-0.00669899,-0.0765414,-0.0683251,-0.0886656 },
        { 0.4207632, 0.3202388, 0.18757344, 0.0830685, 0.0231448,-0.0035101 },
        { 0.2104190, 0.1852464, 0.14351891, 0.0975608, 0.0575074, 0.0282319 },
        { 0.0420838, 0.0410769, 0.03913204, 0.0363792, 0.0329951, 0.0291842 } };
  
      for ( int nl = 1; nl <= 6; ++nl ){
        THEN( "That number of s vector entries are returned, as well as the "
              "energy of the next bragg peak" ){
          std::vector<double> s(nl,0.0), sThisEp(nl);
          for (size_t i = 0; i < epVec.size(); ++i){
            copy(sVariousEp[i].begin(), sVariousEp[i].begin()+nl,sThisEp.begin());
            eNext = computeCrossSections( epVec[i], vec1, vec2, emax, scon, recon,  
                                        s );
            REQUIRE( eNext == Approx(nextBraggEdge[i]).epsilon(1e-6) );
            REQUIRE( ranges::equal( s, sThisEp, equal_1e5 ) );
          }
        } // THEN
      } 
    } // WHEN

    WHEN( "choose beryllium oxide" ){
      int lat = 3, numAtoms = 2;
      double temp  = 650.0, 
             emax  = 8.0,
             recon = 5.1803120897E-20,
             eNext;
      double scon = prepareBraggEdges(lat, temp, emax, numAtoms, vec1, vec2);
  
      using std::copy;

      std::vector<double> epVec {5e-3,1e-2,5e-2,1e-1,5e-1,1,5};
      std::vector<double> nextBraggEdge { 0.007999063, 0.011263090, 0.050768090, 
        0.100972500, 0.501767000, 1.003862000, 5.001879000 };
          
      std::vector<std::vector<double>> sVariousEp {
        { 4.468878,-3.1164850, 1.2491054, 0.0939664,-0.4165726, 0.0251270 },
        { 2.821455,-0.0140015,-0.9892607,-0.2559700, 0.3532396, 0.5913847 },
        { 3.869033, 0.1763646,-0.1249929,-0.3130906,-0.1266984,-0.3194389 },
        { 2.812010, 0.6460590,-0.0658992,-0.0747775,-0.1566545,-0.0726861 },
        { 0.717623, 0.5324232, 0.2982827, 0.1236309, 0.0290206,-0.0106603 },
        { 0.358949, 0.3124967, 0.2368809, 0.1559373, 0.0879284, 0.0402760 },
        { 0.071789, 0.0699317, 0.0663538, 0.0613156, 0.0551662, 0.0483043 } };

  
      for ( int nl = 1; nl <= 6; ++nl ){
        THEN( "That number of s vector entries are returned, as well as the "
              "energy of the next bragg peak" ){
          std::vector<double> s(nl,0.0), sThisEp(nl);
          for (size_t i = 0; i < epVec.size(); ++i){
            copy(sVariousEp[i].begin(), sVariousEp[i].begin()+nl,sThisEp.begin());
            eNext = computeCrossSections( epVec[i], vec1, vec2, emax, scon, recon,  
                                        s );
            REQUIRE( eNext == Approx(nextBraggEdge[i]).epsilon(1e-6) );
            REQUIRE( ranges::equal( s, sThisEp, equal_1e5 ) );
          }
        } // THEN
      } 
    } // WHEN
  } // GIVEN
} // TEST CASE















TEST_CASE( "preparing bragg edges for coherent elastic calculations" ){
  std::vector<double> p(6);
  GIVEN( "graphite is the requested material" ){
    double temp = 296, emax;
    int lat = 1, natom = 1;
    std::vector<double> vec1(5000,0.0), vec2(5000,0.0);


    std::vector<double> vec1_0_100, vec2_0_100, vec1_100_200, vec2_100_200,
                        vec1_200_300, vec2_200_300;
    vec1_0_100 = { 8.794479E15, 8.794479E15, 3.517791E16, 3.517791E16, 
    7.915031E16, 7.915031E16, 8.717302E16, 8.717302E16, 8.717302E16, 9.596750E16, 
    9.596750E16, 9.596750E16, 1.223509E17, 1.223509E17, 1.223509E17, 1.407116E17, 
    1.407116E17, 1.663233E17, 1.663233E17, 1.663233E17, 2.198619E17, 2.198619E17, 
    2.278846E17, 2.278846E17, 2.278846E17, 2.615190E17, 2.615190E17, 2.703135E17, 
    2.703135E17, 2.966969E17, 2.966969E17, 3.070350E17, 3.070350E17, 3.070350E17, 
    3.166012E17, 3.166012E17, 3.406693E17, 3.406693E17, 3.486921E17, 3.486921E17, 
    3.486921E17, 3.574865E17, 3.574865E17, 3.574865E17, 3.838700E17, 3.838700E17, 
    3.838700E17, 4.022307E17, 4.022307E17, 4.037742E17, 4.037742E17, 4.037742E17, 
    4.278424E17, 4.278424E17, 4.278424E17, 4.309294E17, 4.309294E17, 4.813810E17, 
    4.813810E17, 4.894037E17, 4.894037E17, 4.894037E17, 5.181025E17, 5.181025E17, 
    5.181025E17, 5.628466E17, 5.628466E17, 5.685540E17, 5.685540E17, 5.685540E17, 
    5.781203E17, 5.781203E17, 6.102111E17, 6.102111E17, 6.102111E17, 6.190056E17, 
    6.190056E17, 6.190056E17, 6.453890E17, 6.453890E17, 6.453890E17, 6.500197E17, 
    6.500197E17, 6.500197E17, 6.652933E17, 6.652933E17, 6.652933E17, 6.893614E17, 
    6.893614E17, 6.893614E17, 6.924485E17, 6.924485E17, 7.123528E17, 7.123528E17, 
    7.509228E17, 7.509228E17, 7.509228E17, 7.796215E17, 7.796215E17, 7.796215E17 };
    vec2_0_100 = {0.000000E00, 0.000000E00, 2.098463E-8, 2.098463E-8, 
    0.000000E00, 0.000000E00, 1.626950E-9, 1.626950E-9, 1.626950E-9, 9.266133E-9, 
    9.266133E-9, 9.266133E-9, 2.702514E-9, 2.702514E-9, 2.702514E-9, 9.995415E-9, 
    9.995415E-9, 6.814547E-9, 6.814547E-9, 6.814547E-9, 0.000000E00, 0.000000E00, 
    1.886442E-9, 1.886442E-9, 1.886442E-9, 1.387152E-8, 2.774304E-8, 0.000000E00, 
    5.837772E-34,2.562865E-8, 5.125730E-8, 4.701382E-9, 4.701382E-9, 4.701382E-9, 
    6.145991E-9, 6.145991E-9, 0.000000E00, 5.034629E-34,7.213234E-10,7.213234E-10, 
    7.213234E-10,4.257127E-9, 4.257127E-9, 4.257127E-9, 1.352899E-9, 1.352899E-9, 
    1.352899E-9, 4.193765E-8, 2.096882E-8, 1.307116E-9, 1.307116E-9, 1.307116E-9, 
    3.767537E-9, 3.767537E-9, 3.767537E-9, 0.000000E00, 0.000000E00, 0.000000E00, 
    3.970045E-34,1.141440E-9, 1.141440E-9, 1.141440E-9, 3.284509E-9, 3.284509E-9, 
    3.284509E-9, 4.116125E-9, 4.116125E-9, 3.063511E-9, 3.063511E-9, 3.063511E-9, 
    1.613188E-8, 3.226376E-8, 9.670009E-10,9.670009E-10,9.670009E-10,5.737398E-9, 
    5.737398E-9, 5.737398E-9, 1.850387E-9, 1.850387E-9, 1.850387E-9, 9.199318E-10, 
    9.199318E-10,9.199318E-10,9.029481E-10,9.029481E-10,9.029481E-10,5.263708E-9, 
    5.263708E-9, 5.263708E-9, 3.004035E-34,0.000000E00, 0.000000E00, 0.000000E00, 
    1.634200E-9, 1.634200E-9, 1.634200E-9, 2.374227E-9, 2.374227E-9, 2.374227E-9 };
    vec1_100_200 = {7.845572E17, 7.845572E17, 7.845572E17, 7.933517E17, 
    7.933517E17, 7.933517E17, 7.995258E17, 7.995258E17, 7.995258E17, 8.197351E17, 
    8.197351E17, 8.197351E17, 8.243657E17, 8.243657E17, 8.300732E17, 8.300732E17, 
    8.300732E17, 8.637075E17, 8.637075E17, 8.637075E17, 8.794479E17, 8.794479E17, 
    9.115388E17, 9.115388E17, 9.115388E17, 9.252689E17, 9.252689E17, 9.252689E17, 
    9.268124E17, 9.268124E17, 9.268124E17, 9.666210E17, 1.041141E18, 1.061045E18, 
    1.064132E18, 1.125227E18, 1.133249E18, 1.140967E18, 1.151305E18, 1.212400E18, 
    1.215487E18, 1.228140E18, 1.266405E18, 1.347404E18, 1.353578E18, 1.429946E18, 
    1.486267E18, 1.564179E18, 1.573440E18, 1.656287E18, 1.664005E18, 1.723718E18, 
    1.810891E18, 1.909784E18, 1.957615E18, 1.978758E18, 2.087217E18, 2.096478E18, 
    2.107121E18, 2.214504E18, 2.251387E18, 2.368640E18, 2.393480E18, 2.399188E18, 
    2.458900E18, 2.508275E18, 2.541605E18, 2.670273E18, 2.702364E18, 2.710081E18, 
    2.742172E18, 2.763315E18, 2.849411E18, 3.003691E18, 3.018965E18, 3.035944E18, 
    3.174807E18, 3.358091E18, 3.373526E18, 3.517792E18, 3.701076E18, 3.866484E18, 
    3.878365E18, 4.077390E18, 4.082020E18, 4.128003E18, 4.139884E18, 4.158083E18, 
    4.227057E18, 4.256528E18, 4.488577E18, 4.491340E18, 4.569575E18, 4.605220E18, 
    4.652280E18, 4.886015E18, 4.912560E18, 5.000972E18, 5.065620E18, 5.326349E18 };
    vec2_100_200 = { 6.297022E-09, 6.297022E-09, 6.297022E-09, 1.339646E-34, 
    1.33964E-34, 5.35858E-34, 2.32313E-09, 2.32313E-09, 2.32313E-09, 1.21232E-08, 
    1.21232E-08, 1.21232E-08, 2.41267E-08, 1.20633E-08, 4.49638E-09, 4.49638E-09, 
    4.49638E-09, 1.24306E-34, 1.24306E-34, 4.97224E-34, 2.84686E-09, 2.84686E-09, 
    6.88838E-10, 6.88838E-10, 6.88838E-10, 1.08705E-08, 1.08705E-08, 1.08705E-08, 
    1.35671E-09, 1.35671E-09, 1.35671E-09, 1.95659E-09, 2.54375E-08, 5.36451E-09, 
    5.56516E-08, 7.72558E-34, 1.67377E-09, 3.65587E-08, 3.69482E-08, 9.36236E-09, 
    5.34831E-34, 1.53918E-09, 1.54523E-08, 3.05644E-08, 3.19697E-08, 7.62341E-09, 
    2.29103E-08, 7.01180E-09, 2.60367E-08, 1.08858E-09, 3.55223E-08, 1.08708E-08, 
    6.48701E-08, 1.42732E-08, 3.02210E-08, 3.16336E-08, 4.77266E-09, 4.74187E-09, 
    3.37862E-08, 2.18044E-09, 2.56310E-08, 1.42925E-08, 2.06481E-08, 7.56813E-09, 
    1.33323E-08, 2.38384E-08, 1.55647E-08, 8.60645E-09, 5.26869E-10, 2.40743E-08, 
    1.53062E-09, 1.20400E-08, 2.17052E-08, 8.70192E-10, 8.61914E-10, 4.87628E-08, 
    2.73442E-08, 6.95854E-10, 3.07499E-08, 3.92303E-08, 2.63412E-08, 2.57924E-10, 
    2.42334E-08, 3.64728E-09, 7.72361E-09, 4.42686E-10, 2.51825E-35, 1.00049E-08, 
    8.33930E-10, 2.72397E-08, 1.07905E-09, 1.23403E-08, 3.71767E-09, 1.83415E-09, 
    2.38537E-08, 8.34468E-34, 5.75221E-09, 1.00928E-08, 2.10682E-08, 6.73928E-10 };

    vec1_200_300 = { 5.327139E18, 5.335610E18, 5.354133E18, 5.414312E18, 
    5.496550E18, 5.785529E18, 5.798936E18, 5.845242E18, 5.945068E18, 6.249159E18, 
    6.250236E18, 6.281107E18, 6.293760E18, 6.411175E18, 6.743194E18, 6.759868E18, 
    6.894872E18, 7.243079E18, 7.243564E18, 7.396157E18, 7.766117E18, 7.767984E18, 
    7.773368E18, 7.789270E18, 7.790508E18, 7.805944E18, 7.915031E18, 8.326953E18, 
    8.327420E18, 8.428485E18, 8.451495E18, 8.876537E18, 8.886713E18, 8.896298E18, 
    8.905235E18, 8.923758E18, 9.005547E18, 9.457439E18, 9.549404E18, 9.575483E18, 
    9.577188E18, 1.006148E19, 1.009436E19, 1.016642E19, 1.068109E19, 1.068265E19, 
    1.077324E19, 1.131478E19, 1.131754E19, 1.139765E19, 1.197358E19, 1.197481E19, 
    1.199025E19, 1.203964E19, 1.264985E19, 1.269363E19, 1.269923E19, 1.335106E19, 
    1.336958E19, 1.337640E19, 1.405216E19, 1.406573E19, 1.407117E19, 1.478352E19, 
    1.554352E19, 1.555573E19, 1.556808E19, 1.560063E19, 1.626099E19, 1.708591E19, 
    1.711328E19, 1.780882E19, 1.870184E19, 1.870598E19, 1.880963E19, 1.888372E19, 
    1.895548E19, 1.895781E19, 1.942700E19, 2.040281E19, 2.042586E19, 2.043480E19, 
    2.043742E19, 2.053247E19, 2.061117E19, 2.111554E19, 2.217406E19, 2.217826E19, 
    2.220710E19, 2.224772E19, 2.231629E19, 2.233489E19, 2.287444E19, 2.403246E19, 
    2.404001E19, 2.408713E19, 2.412896E19, 2.470369E19, 2.594474E19, 2.595960E19 };
    vec2_200_300 = { 1.79635E-09, 4.77775E-09, 3.14648E-09, 4.04778E-09, 
    1.93835E-08, 5.23574E-10, 3.27535E-09, 9.00499E-09, 1.47686E-08, 2.17105E-09, 
    1.14758E-09, 5.97140E-10, 4.33414E-09, 1.20139E-08, 1.75607E-09, 2.88790E-09, 
    1.46760E-08, 7.98099E-11, 4.45305E-09, 9.40212E-09, 6.06022E-11, 6.05429E-11, 
    9.64777E-10, 8.98070E-11, 4.75194E-10, 2.72496E-09, 8.26345E-09, 1.59554E-34, 
    2.78215E-09, 3.62467E-10, 6.43188E-09, 3.40224E-11, 1.01532E-10, 5.00307E-35, 
    1.00567E-10, 1.27815E-09, 5.89880E-09, 3.43224E-10, 6.19969E-10, 3.56326E-11, 
    5.32238E-09, 9.24668E-11, 3.85898E-10, 3.98730E-09, 1.08237E-10, 4.33485E-10, 
    3.04279E-09, 1.57170E-10, 3.83569E-10, 2.37407E-09, 7.05390E-12, 4.72475E-35, 
    2.08084E-10, 1.75703E-09, 6.49758E-11, 3.93617E-11, 1.53321E-09, 1.06384E-11, 
    5.27048E-12, 1.16342E-09, 4.00668E-11, 1.24370E-12, 7.84610E-10, 6.58829E-10, 
    9.59600E-12, 2.98028E-12, 1.05992E-11, 3.72277E-10, 3.57909E-10, 2.41781E-11, 
    2.17721E-10, 1.64309E-10, 7.68023E-13, 2.56775E-11, 9.86660E-12, 7.92668E-12, 
    5.65656E-13, 5.90950E-11, 9.37659E-11, 6.72805E-13, 2.10614E-12, 1.10411E-13, 
    1.12497E-11, 2.99441E-12, 3.19945E-11, 4.63431E-11, 3.64133E-37, 9.49684E-13, 
    3.70005E-13, 1.50062E-12, 3.43035E-13, 1.64378E-11, 1.87597E-11, 1.55808E-13, 
    5.28610E-13, 1.12779E-12, 7.57311E-12, 9.76623E-12, 7.78169E-15, 1.25947E-13 };


    WHEN( "smal value for emax" ){
      emax = 0.05;
      THEN ( "Not many bragg edges, rest of vec1 and vec2 are zero" ){
        prepareBraggEdges( lat, temp, emax, natom, vec1, vec2 );
        REQUIRE( vec1.size() == 132 );
        REQUIRE( vec2.size() == 132 );

        for ( size_t i = 0; i < 100; ++i ){ 
          REQUIRE( vec1_0_100[i] == Approx( vec1[i] ).epsilon(1e-6) );
          REQUIRE( vec2_0_100[i] == Approx( vec2[i] ).epsilon(1e-6) );
        }

        for ( size_t i = 100; i < vec1.size()-1; ++i ){ 
          REQUIRE( vec1_100_200[i-100] == Approx( vec1[i] ).epsilon(1e-6) );
          REQUIRE( vec2_100_200[i-100] == Approx( vec2[i] ).epsilon(1e-6) );
        }

      } // THEN 
    } // WHEN

  
    WHEN( "Medium value for emax" ){
      emax = 1.2;
      THEN ( "Healthy amount of bragg edges, rest of vec1 and vec2 are zero" ){
        prepareBraggEdges( lat, temp, emax, natom, vec1, vec2 );
        REQUIRE( vec1.size() == 294 );
        REQUIRE( vec2.size() == 294 );

        for ( size_t i = 0; i < 100; ++i ){ 
          REQUIRE( vec1_0_100[i] == Approx( vec1[i] ).epsilon(1e-6) );
          REQUIRE( vec2_0_100[i] == Approx( vec2[i] ).epsilon(1e-6) );
        }

        for ( size_t i = 100; i < 200; ++i ){ 
          REQUIRE( vec1_100_200[i-100] == Approx( vec1[i] ).epsilon(1e-6) );
          REQUIRE( vec2_100_200[i-100] == Approx( vec2[i] ).epsilon(1e-6) );
        }

        for ( size_t i = 200; i < vec1.size()-1; ++i ){ 
          REQUIRE( vec1_200_300[i-200] == Approx( vec1[i] ).epsilon(1e-6) );
          REQUIRE( vec2_200_300[i-200] == Approx( vec2[i] ).epsilon(1e-6) );
        }

      } // THEN 
    } // WHEN
    WHEN( "large value for emax" ){
      emax = 5.5;
      THEN ( "Respectable # of bragg edges, rest of vec1 and vec2 are zero" ){
        prepareBraggEdges( lat, temp, emax, natom, vec1, vec2 );
        REQUIRE( vec1.size() == 435 );
        REQUIRE( vec2.size() == 435 );

        for ( size_t i = 0; i < 100; ++i ){ 
          REQUIRE( vec1_0_100[i] == Approx( vec1[i] ).epsilon(1e-6) );
          REQUIRE( vec2_0_100[i] == Approx( vec2[i] ).epsilon(1e-6) );
        }

        for ( size_t i = 100; i < 200; ++i ){ 
          REQUIRE( vec1_100_200[i-100] == Approx( vec1[i] ).epsilon(1e-6) );
          REQUIRE( vec2_100_200[i-100] == Approx( vec2[i] ).epsilon(1e-6) );
        }

       for ( int i = 200; i < 300; ++i ){ 
          REQUIRE( vec1_200_300[i-200] == Approx( vec1[i] ).epsilon(1e-6) );
          REQUIRE( vec2_200_300[i-200] == Approx( vec2[i] ).epsilon(1e-6) );
        }

      } // THEN 
    } // WHEN
  } // GIVEN






  GIVEN( "beryllium is the requested material" ){
    double temp = 500.0, emax;
    int lat = 2, natom = 1;
    std::vector<double> vec1(5000,0.0), vec2(5000,0.0);


    std::vector<double> vec1_0_100, vec2_0_100, vec1_100_200, vec2_100_200,
                        vec1_200_300, vec2_200_300;
    vec1_0_100 = {3.074805E16, 3.074805E16, 1.007623E17, 1.007623E17, 
    1.007623E17, 1.229922E17, 1.229922E17, 1.315103E17, 1.315103E17, 1.315103E17, 
    2.237545E17, 2.237545E17, 2.237545E17, 2.767325E17, 2.767325E17, 3.022869E17, 
    3.022869E17, 3.330349E17, 3.330349E17, 3.774948E17, 3.774948E17, 3.774948E17, 
    4.030492E17, 4.030492E17, 4.030492E17, 4.252791E17, 4.252791E17, 4.337972E17, 
    4.337972E17, 4.337972E17, 4.919689E17, 4.919689E17, 5.260414E17, 5.260414E17, 
    5.260414E17, 5.790194E17, 5.790194E17, 5.927312E17, 5.927312E17, 5.927312E17, 
    6.797816E17, 6.797816E17, 6.797816E17, 7.053360E17, 7.053360E17, 7.053360E17, 
    7.360841E17, 7.360841E17, 7.360841E17, 7.687014E17, 7.687014E17, 7.942557E17, 
    7.942557E17, 8.283282E17, 8.283282E17, 8.283282E17, 8.694636E17, 8.694636E17, 
    8.694636E17, 8.950180E17, 8.950180E17, 8.950180E17, 9.068606E17, 9.068606E17, 
    9.068606E17, 9.376086E17, 9.376086E17, 9.376086E17, 1.029853E18, 1.106930E18, 
    1.171751E18, 1.207692E18, 1.309910E18, 1.398829E18, 1.474037E18, 1.506655E18, 
    1.586642E18, 1.607417E18, 1.701116E18, 1.735189E18, 1.888929E18, 1.909704E18, 
    1.967875E18, 2.068638E18, 2.191216E18, 2.211991E18, 2.270162E18, 2.370925E18, 
    2.490592E18, 2.642049E18, 2.719127E18, 2.874736E18, 2.893642E18, 3.074805E18, 
    3.246623E18, 3.277785E18, 3.287759E18, 3.397453E18, 3.477855E18, 3.658190E18};
    vec2_0_100   = { 0.000000E00, 0.000000E00, 2.857665E-9, 2.857665E-9, 
    2.857665E-9, 5.063031E-9, 5.063031E-9, 1.456839E-8, 1.456839E-8, 1.456839E-8, 
    3.405052E-9, 3.405052E-9, 3.405052E-9, 0.000000E00, 0.000000E00, 5.430392E-9, 
    1.086078E-8, 0.000000E00, 0.000000E00, 6.777581E-9, 6.777581E-9, 6.777581E-9, 
    1.066503E-9, 1.066503E-9, 1.066503E-9, 1.625860E-8, 8.129302E-9, 5.987270E-9, 
    5.987270E-9, 5.987270E-9, 1.771490E-9, 1.771490E-9, 1.657603E-9, 1.657603E-9, 
    1.657603E-9, 0.000000E00, 0.000000E00, 1.463993E-9, 1.463993E-9, 1.463993E-9, 
    3.769870E-9, 3.769870E-9, 3.769870E-9, 1.203522E-9, 1.203522E-9, 1.203522E-9, 
    6.861505E-9, 6.861505E-9, 6.861505E-9, 0.000000E00, 0.000000E00, 4.162634E-9, 
    8.325269E-9, 1.971969E-9, 1.971969E-9, 1.971969E-9, 2.774483E-9, 2.774483E-9, 
    2.774483E-9,8.892684E-10,8.892684E-10,8.892684E-10, 1.746756E-9, 1.746756E-9, 
    1.746756E-9, 0.000000E00, 1.404663E-8, 0.000000E00, 8.731435E-9, 1.302790E-9, 
    8.795036E-9, 5.084491E-9, 1.570423E-8, 1.291964E-8, 7.123059E-9, 1.132874E-9, 
    6.156946E-9, 6.359091E-9, 3.548504E-9, 3.957397E-9, 2.105949E-9, 6.663487E-9, 
    4.005657E-9, 6.575435E-9, 2.918932E-9, 6.523584E-9, 3.448611E-9, 3.057076E-9, 
    4.416373E-9, 2.477778E-9, 5.402565E-9,8.769307E-10, 1.607555E-9, 6.233046E-9,
    2.879079E-10,2.780258E-10,1.177223E-9, 1.438005E-9, 1.273365E-9, 0.000000E00 };
    vec1_100_200 = {3.720515E18, 3.928274E18, 3.981666E18, 4.123564E18, 
    4.332778E18, 4.384715E18, 4.396144E18, 4.427720E18, 4.687002E18, 4.830769E18, 
    5.091506E18, 5.133056E18, 5.134097E18, 5.196421E18, 5.498708E18, 5.593863E18, 
    5.599470E18, 5.896149E18, 5.901757E18, 6.026619E18, 6.328905E18, 6.429668E18, 
    6.781821E18, 6.808618E18, 6.918312E18, 7.312429E18, 7.321361E18, 7.688682E18, 
    7.711326E18, 7.715478E18, 7.825173E18, 7.871502E18, 8.274551E18, 8.718949E18, 
    8.760498E18, 8.778362E18, 8.886188E18, 9.365072E18, 9.437369E18, 9.483698E18, 
    9.793048E18, 9.962370E18, 1.047159E19, 1.049838E19, 1.086923E19, 1.110005E19, 
    1.165996E19, 1.177007E19, 1.180538E19, 1.185566E19, 1.200691E19, 1.229922E19, 
    1.291980E19, 1.295096E19, 1.296094E19, 1.315103E19, 1.320608E19, 1.355989E19, 
    1.425942E19, 1.426523E19, 1.432029E19, 1.446675E19, 1.488206E19, 1.563726E19, 
    1.567590E19, 1.578892E19, 1.626572E19, 1.708076E19, 1.712190E19, 1.715182E19, 
    1.717258E19, 1.771088E19, 1.859781E19, 1.861774E19, 1.921753E19, 2.019545E19, 
    2.022994E19, 2.078568E19, 2.183113E19, 2.184630E19, 2.200917E19, 2.204366E19, 
    2.209559E19, 2.215461E19, 2.237545E19, 2.239788E19, 2.241533E19, 2.355031E19, 
    2.363612E19, 2.402753E19, 2.410647E19, 2.531562E19, 2.535781E19, 2.541638E19, 
    2.560982E19, 2.566632E19, 2.571867E19, 2.585911E19, 2.716902E19, 2.723447E19 };
    vec2_100_200 = { 2.771224E-09, 4.056048E-10, 2.128007E-09, 8.063332E-10, 
    2.964131E-10, 8.237195E-11, 7.946766E-10, 7.770736E-10, 4.796811E-10, 
    9.009425E-10, 1.530877E-10, 3.690742E-11, 1.388272E-10, 7.622018E-10, 
    7.371049E-11, 5.652275E-11, 5.539253E-10, 8.218794E-11, 1.242771E-10, 
    2.148807E-10, 9.036160E-11, 2.492934E-10, 3.863648E-11, 1.268326E-10, 
    1.134038E-10, 0.000000E00, 1.141030E-10, 7.631996E-12, 9.940848E-12, 
    1.333637E-11, 4.295375E-12, 7.877221E-11, 4.601739E-11, 5.212715E-12, 
    1.686038E-12, 1.025461E-11, 2.373603E-11, 3.304939E-12, 4.380227E-12, 
    1.019992E-11, 6.949444E-12, 1.037046E-11, 1.029450E-12, 4.393329E-12, 
    1.554998E-12, 4.001470E-12, 8.587273E-13, 4.683730E-13, 1.147388E-13, 
    6.691252E-13, 8.573953E-13, 1.293856E-12, 1.243738E-14, 1.205346E-14, 
    2.234719E-13, 7.369851E-14, 2.712956E-13, 4.782004E-13, 1.295589E-14, 
    4.115180E-14, 3.564478E-14, 8.781616E-14, 1.338788E-13, 2.446460E-15, 
    8.406169E-15, 4.504527E-14, 2.991807E-14, 1.930601E-16, 1.477330E-15, 
    2.697896E-16, 1.192060E-14, 8.858012E-15, 2.979180E-16, 2.963182E-15, 
    2.006583E-15, 3.426932E-17, 5.707316E-16, 5.253637E-16, 5.169554E-18, 
    2.979171E-17, 9.885421E-18, 1.224696E-17, 1.794045E-17, 3.447411E-17, 
    5.487573E-18, 3.399470E-18, 1.114770E-16, 5.150730E-18, 1.192204E-17, 
    4.442415E-18, 2.083696E-17, 2.685393E-19, 6.204995E-19, 1.257043E-18, 
    3.166091E-19, 4.012800E-19, 6.480621E-19, 4.327015E-18, 1.378948E-19, 
    1.135265E-19};

    WHEN( "small value for emax" ){
      emax = 0.8;
      THEN ( "Not many bragg edges, rest of vec1 and vec2 are zero" ){
        prepareBraggEdges( lat, temp, emax, natom, vec1, vec2 );
        REQUIRE( vec1.size() == 164 );
        REQUIRE( vec2.size() == 164 );

        for ( size_t i = 0; i < 100; ++i ){ 
          REQUIRE( vec1_0_100[i] == Approx( vec1[i] ).epsilon(1e-6) );
          REQUIRE( vec2_0_100[i] == Approx( vec2[i] ).epsilon(1e-6) );
        }

        for ( size_t i = 100; i < vec1.size()-1; ++i ){ 
          REQUIRE( vec1_100_200[i-100] == Approx( vec1[i] ).epsilon(1e-6) );
          REQUIRE( vec2_100_200[i-100] == Approx( vec2[i] ).epsilon(1e-6) );
        }

      } // THEN 
    } // WHEN

  
    WHEN( "Medium value for emax" ){
      emax = 2.5;
      THEN ( "Healthy amount of bragg edges, rest of vec1 and vec2 are zero" ){
        prepareBraggEdges( lat, temp, emax, natom, vec1, vec2 );
        REQUIRE( vec1.size() == 250 );
        REQUIRE( vec2.size() == 250 );

        for ( size_t i = 0; i < 100; ++i ){ 
          REQUIRE( vec1_0_100[i] == Approx( vec1[i] ).epsilon(1e-6) );
          REQUIRE( vec2_0_100[i] == Approx( vec2[i] ).epsilon(1e-6) );
        }

        for ( size_t i = 100; i < 200; ++i ){ 
          REQUIRE( vec1_100_200[i-100] == Approx( vec1[i] ).epsilon(1e-6) );
          REQUIRE( vec2_100_200[i-100] == Approx( vec2[i] ).epsilon(1e-6) );
        }

      } // THEN 
    } // WHEN
    WHEN( "large value for emax" ){
      emax = 5.0;
      THEN ( "Respectable # of bragg edges, rest of vec1 and vec2 are zero" ){
        prepareBraggEdges( lat, temp, emax, natom, vec1, vec2 );
        REQUIRE( vec1.size() == 296 );
        REQUIRE( vec2.size() == 296 );
  
        for ( size_t i = 0; i < 100; ++i ){ 
          REQUIRE( vec1_0_100[i] == Approx( vec1[i] ).epsilon(1e-6) );
          REQUIRE( vec2_0_100[i] == Approx( vec2[i] ).epsilon(1e-6) );
        }

        for ( size_t i = 100; i < 200; ++i ){ 
          REQUIRE( vec1_100_200[i-100] == Approx( vec1[i] ).epsilon(1e-6) );
          REQUIRE( vec2_100_200[i-100] == Approx( vec2[i] ).epsilon(1e-6) );
        }

      } // THEN 
    } // WHEN
  } // GIVEN
} // TEST CASE
